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The  objective  of  this  research  was  to  systematically  explore  the  attributes  of  the 
wave  form  Lg  and  related  wave  forms  which  are  trapped  in  the  crust.  The  different 
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was  achieved  through  waveform  data  analyses.  The  Kennett  conjecture  that  Lg  waves 
can  be  adequately  described  by  multiple  S  energy  trapped  in  the  crust  has  been 
adopted.  It  has  been  shown  that  Lg  magnitudes  are  much  more  stable  than  the  chaotic 
behavior  of  the  ray  would  seem  to  imply.  As  a  result  of  recognizing  that  a  large 
number  of  caustic  are  generated  as  a  consequence  of  the  chaotic  ray  beharior,  it  was 
decided  that  Maslov  theory  was  probably  the  only  suitable  technique  to  deal  with  the 
singularity  problem.  As  a  result.  It  has  been  shown  that  Maslov  synthetics  can  be 
used  to  justify  the  use  of  ray  density  plots  as  energy  plots,  which  can  be  used  in 
turn  for  the  calibration  of  Lg  magnitudes.  The  3-component  array  and  single  station 
data  from  the  NPE  1993  detonation  have  been  analyzed  in  order  to  examine  the  effects 
boundary  undulations  of  the  crustal  wave  guide  may  have  on  propagation. 
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Our  project  aimed  for  a  systematic  exploration  of  the  attributes  of  Lg  and 
related  trapped  waves  in  the  crust.  We  studied  the  different  theoretical 
approximations  used  to  characterize  or  fully  compute  regional  seismic  waves, 
and  analyzed  waveform  data  to  develop  an  empirical  classification  of  the 
behaviour  of  Lg. 

From  the  beginning  we  adopted  Kennett's  conjecture  that  Lg  waves  could 
be  adequately  described  by  multiple  S  energy  trapped  in  the  crust.  We  de¬ 
veloped  software  to  generate  synthetics  by  ray  summation  and  validated  this 
against  a  mode-summation  code  for  horizontally  layered  models  [1-4,7,14]. 
We  then  .set  out  to  perturb  the  depth  to  the  Moho  to  obtain  more  realistic, 
complicated  seismograms  [12].  When  we  discovered  that  the  chaotic  nature 
of  the  kinematics  of  the  trapped  rays  is  a  major  characteristic  of  Lg  waves 
[5.8],  we  judged  this  discovery  important  enough  to  warrant  extra  attention. 
In  a  first  application,  we  showed  that  Lg  magnitudes  are  much  more  stable 
measurements  than  the  chaotic  behaviour  of  the  ray  would  led  one  to  expect 
[5.11] 

One  consequence  of  the  chaotic  ray  behavior  is  the  large  numbei  of  caustics 
that  are  generated.  These,  in  turn,  lead  to  singularities  that  render  the  more 
common  ray  summation  techniques  useless.  We  decided  that  Maslov  theory 
WcLS  probably  the  only  suitable  technique  to  deal  with  the  singularity  problem 
[6],  and  set  out  to  compute  the  correct  Maslov  phase  and  amplitude  using 
the  method  of  stationary  phase.  We  called  in  the  help  of  Chris  Chapman  at 
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Schlumberger  (UK),  and  Henk  Keers  spent  the  summer  of  1995  working  with 
him  in  Cambridge  to  solve  the  problems  posed  by  pseudo-caustics  in  calcu¬ 
lating  Maslov  synthetics  (pseudo  caustics  are  essentially  plane  wave  arrivals 
that  often  accompany  regular  caustic  singularities)  [6,10,11]. 

An  algorithm  is  ready  that  efficiently  computes  2D  Maslov  integrals.  This 
algorithm  is  essential  for  the  fast  computation  of  synthetics  in  3D  heteroge¬ 
neous  media.  The  latter  development  is  still  in  progress,  and  will  form  the 
concluding  chapter  of  Keers’  PhD  thesis.  Very  complete  seismograms  can  be 
generated  in  crustal  models  of  arbitrary  complexity  through  the  generation 
of  multiples,  for  w’hich  we  developed  an  efficient  tree  algorithm  [8]. 

In  a  parallel  observational  study,  we  made  sure  that  our  theoretical  foun¬ 
dations  where  always  strongly  based  on  the  reality  of  regional  seismograms. 
This  part  of  the  project  was  carried  by  Kristin  Vogfjord.  She  analyzed  3- 
component  array-  and  single-station  data,  recording  the  NPE  explosion  in 
1993.  on  a  profile  extending  from  the  NTS,  westward  across  the  Sierra  Moun¬ 
tains  [18.  21.  23].  PmP  and  SmS  are  coherent  into  the  Sierras,  but  disappear 
due  to  the  crustal  gradient  at  around  230  km  distance.  They  arrive  off- 
azimuth  indicating  a  dip  on  the  Moho  at  the  reflection.  The  second  multiple 
•Molio  reflections  are  not  so  clearly  observed.  Due  to  the  Moho  undulations 
tlie  rays  in  the  second  Moho  multiple  start  to  diverge  and  are  not  clearly 
traced  into  the  Sierras.  Only  one  of  the  arrays  has  coherent  arrivals  corre¬ 
sponding.  in  time  and  phase  velocity,  to  2PmP  and  2SmS.  This  confirmed 
our  theoretical  findings  that,  while  multiple  S  is  important  in  building  up 
the  Lg  wavetrain.  internal  multiples  and  crustal  heterogeneity  are  important 
and  must  be  coinsidered  when  explaining  Lg  and  its  characteristics. 

The  observation  that  S  is  important  even  in  explosively  generated  waves  re¬ 
quires  a  theoretical  understanding.  In  the  literature,  there  is  no  consensus  on 
this  topic.  X’ogfjord  established  that  explosions  can  generate  sufficient  S*  en¬ 
ergy  to  allow  for  the  validity  of  Kennett’s  hypothesis  [13,15,16,24].  Vogfjord 
examined  effects  of  two  parameters,  source  depth  and  Moho  undulations,  on 
the  development  and  character  of  the  crustal  wave  trains  (P,  Lp).  If  source 
depth  is  within  1  km  from  the  surface,  and  the  velocity  structure  is  such 
that  the  peak  in  the  S’  radiation  pattern  is  trapped  in  the  crust,  then  S* 
contributes  significantly  to  Lg.  However  if  the  source  radiates  shear  waves, 
includes  spall,  or  tectonic  release  this  contribution  may  be  drowned  in  the 
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S  waves  emanating  from  these  other  terms.  Explosion  sources  that  do  not 
include  any  significant  amount  of  these  contributing  factors,  on  the  other 
hand,  may  lead  to  S*  dominated  Lg  wave  trains.  Such  sources  are  possibly 
representative  of  decoupled  nuclear  explosions,  as  well  as  some  singe-shot  un¬ 
derground  mining  explosions.  Since  decoupled  nuclear  explosions  represent 
a  plausible  threat  to  a  comprehensive  test  ban,  observations  of  S*  signature 
in  Lg  could  possibly  be  an  important  discriminant  of  such  evasive  testing 
[17,19,24]. 

The  importance  of  the  depth  parameter  for  discrimination  led  to  a  second 
stud}’  on  2Pg  in  the  crustal  P  wave  train  and  2Sg  in  the  Lg  wave  train,  using 
data  from  the  Iceland  SIL  network  [16,  20,  22].  Their  relative  arrival-time  is 
strongly  depth  dependent  and  thus  provides  a  fixed  constraint  on  the  depth 
of  shallow  sources.  The  results  have  been  reported  in  [16],  and  a  paper  is  in 
preparation  [22]. 

The  Maslov  synthetics  can  be  used  to  justify  the  use  of  ray  density  plots 
as  energy  plots,  which  can  be  used  in  turn  for  the  calibration  of  Lg  magni¬ 
tudes.  In  order  to  examine  the  effects  that  undulations  on  the  boundaries 
of  the  crustal  wave-guide  have  on  the  propagation  of  crustal  wave  trains, 
w(‘  analyzed  the  3-component  array-  and  single-station  data  from  the  NPE 
explosion  in  1993.  A  paper  detailing  the  results  is  in  preparation  [23]. 
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The  objective  of  our  research  is  to  find  ways  to  compute  synthetics  of  Lg-waves  in  com¬ 
plicated  heterogeneous  media.  Lg-waves  arc  the  most  prominent  phases  on  regional  seismo¬ 
grams.  As  such  they  are  important  in  studying  seismic  sources  and  the  structure  of  the  crust 
and  upper  mantle. 

In  our  study  we  look  at  Lg-waves  from  both  a  mode  and  a  ray  point  of  view  (Kennett, 
1986).  The  validation  of  the  ray  methods,  by  comparing  them  to  the  mode  sum  is  the  pri¬ 
mary  reason  for  this  study.  Both  methods  have  been  among  the  most  successful  in  computing 
synthetics  in  heterogeneous  media.  However  the  ray  method  tends  to  be  less  difficult  from 
a  computational  point  of  view. 

To  gain  physical  understanding  we  start  by  considering  the  propagation  of  Lg-waves  in 
a  simple  layer  over  a  halfspace.  A  comparison  of  the  seismograms,  computed  using  the 
two  different  methods  is  presented  and  an  interpretation  of  the  Lg-waves  in  terms  of  modes 
and  rays  18  given.  Subsequently  we  discuss  the  effects  of  an  undulating  Moho  on  the  wave 
propagation;  in  doing  so  we  emphasize  ray  theory,  because  of  its  computational  advantages 
over  mode  theory.  First  we  discuss  one-point  ray  tracing  and  point  out  the  limitations  of 
two-point  ray  tracing.  Because  of  the  limitations  of  classical  ray  theory  we  also  study  the 
application  of  graph  theory  to  this  kind  of  models. 

RESEARCH  ACCOMPLISHED 
Lg-waves  in  simple  layered  media 

In  order  to  understand  the  ray-mode  duality  of  Lg-waves  in  layered  media  we  take  a  simple 
layer  over  a  haHspace  model,  see  table  1.  The  source  that  we  use  is  a  strike-slip  point  source 
at  8  km  depth.  AU  seismograms  are  calculated  for  points  on  a  line  through  the  source,’ 
perpendicular  to  the  strike  direction.  For  this  preliminary  investigation  we  only  study  SH 
waves.  The  calculated  synthetics  are  all  convolved  with  the  WWSSN  short  period  instru¬ 
ment  response,  which  acts  as  a  passband  with  comer  frequencies  of  IHz  and  3Hz. 

Modes  -  Using  the  locked  mode  approximation  of  normal  mode  theory,  synthetics  have  been 
calculated.  The  epicentral  distance  was  500  km.  Figure  1  shows  seismograms  in  which  the 
first  ten,  twenty,  forty  and  hundred  modes  are  used.  After  comparing  the  seismograms  we 
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conclude  th&t  ten  modes  are  cleaxly  not  enough  to  model  the  whole  Lg  wavetrain.  If  twenty 
modes  are  used  then  the  first  half  of  the  seismogram  is  reasonably  well  modelled,  when  com¬ 
pared  with  the  seismogram  in  which  hundred  modes  have  been  used.  If  forty  modes  are  used 
then  the  result  is  quite  satisfactory)  except  for  the  fact  that  the  last  part  of  the  seismogram 
is  not  well  modelled.  It  appears  that  we  need  around  sixty  additional  modes  in  order  to  get 
one  more  arrival. 

The  latest  arriving  energy  of  each  mode  corresponds  with  the  Airy  phase.  Even  though 
the  Airy  phases  are  associated  with  a  concentration  of  energy,  the}’  are  not  visible  on  the 
seismograms  as  separate  arrivals,  as  suggested  for  example  by  Kovach  and  Anderson  (1964). 
Rather,  the  formation  of  the  different  multiples  is  a  consequence  of  a  more  complex  super¬ 
position  of  the  various  modes. 

Rays  -  In  a  simple  layer  over  halfspace  model  it  is  easy  to  calculate  all  the  ray  paths  from  the 
source  to  the  receiver,  as  well  as  the  traveltimes  and  amplitudes  of  ray.  The  amplitudes 
depend  on  the  radiation  pattern  of  the  source,  the  pathlength  of  the  ray  and  the  reflection 
coefficient  of  the  Moho.  Using  these  quantities  ray  synthetics  can  readily  be  computed.  An 
example  of  such  synthetics  are  shown  in  figures  2  and  3.  The  only  rays  that  contribute 
significantly  to  the  seismogram  are  the  super-critically  reflected  rays.  In  this  case  (epicentral 
distance  of  500  km)  there  are  twenty-one  super-critical  reflected  rays.  The  twenty-second 
multiple  has  an  amplitude  that  is  a  factor  of  0.7  x  10"^  smaller  than  the  twenty-first  multiple, 
due  to  leakage  of  energy  into  the  mantle.  All  super-critical  reflected  rays  are  clearly  visible 
as  separate  arrivals,  providing  a  simple  interpretation  of  the  seismogram  in  terms  of  rays. 
Because  the  reflection  coefficient  is  complex  the  multiples  have  different  phases.  It  should 
be  noted  that  head-wave  energy  is  not  modelled  with  this  method. 

Mode-Ray  duality  -  In  figure  2  the  mode  and  ray  seismograms  for  different  epicentral 
distances  are  presented.  The  overall  agreement  between  the  two  methods  is  very  good.  In 
figure  3  blow-ups  of  the  seismograms  for  an  epicentral  distance  of  five  hundred  kilometers 
are  given.  It  appears  that  there  are  two  major  differences;  the  first  difference  is  the  pres¬ 
ence  of  energy  in  the  mode  synthetic  between  the  multiples.  This  is  for  instance  visible  in 
figure  3c,  between  166s  and  1688,  and  between  169s  amd  ITls.  This  energy  is  not  present  in 
the  ray  synthetics.  The  second  difference  is  the  last  multiple.  This  multiple  is  not  present 
in  the  mode  summation.  Currently  we  believe  that  the  energy  in  between  the  multiples  is 
associated  with  (multiple)  head-waves.  Head-wave  and  multiple  head-wave  arrivals  are  not 
accounted  for.  The  second  discrepancy  can  be  explained  either  by  the  incompleteness  of  the 
mode  sum,  such  that  more  than  100  modes  need  to  be  taken  into  account,  or  by  the  failure 
of  ray  theory  for  rays  reflected  near  critical  angles. 

It  should  be  noted  that  in  more  complicated  layered  media  not  all  ray  arrivals  are  visible 
as  separate  arrivals,  thus  making  the  interpretation  of  the  seismogram  in  terms  of  rays  more 
complicated. 

Lg-waves  in  simple  heterogeneous  media 

As  a  next  step  we  study  the  behavior  of  Lg-waves  in  a  layer  over  a  halfspace  model  with 
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an  undulating  Moho.  Rrom  a  modal  point  of  view,  in  the  high  frequency  approximation  the 
waveforms  should  be  identical  to  the  model  with  a  flat  Moho  and  a  crustal  equal 

to  the  average  crustal  thickness  between  source  and  receiver.  If  the  high  frequency  approxi¬ 
mation  is  not  made,  then  the  modes  couple  and  the  computation  of  synthetics  becomes  more 
difficult  (Maupin,  1990). 

Rays  -  Prom  a  ray  point  of  view  the  first  thing  to  do  is  one-point  ray  tracing.  The  situation 
is  sketched  in  figure  4.  The  Moho  is  described  by  a  given  function  /.  The  thirlm^sB  of  the 
crust  at  the  horizontal  coordinate  x  is  given  by  h-f/(x).  li,  after  n  reflections  the  horizontal 
coordinate  of  the  ray  is  x„  then  after  n  -|- 1  reflections  the  horizontal  coordinate  of  the  ray  is 
given  by 

2n+l  =  V'n  +  [A  + /(V’n)]/ tan  tfn+l  (1) 

where  V'n+:  and  6„+i  are  given  by 

V’n  =  x„-f  [h  +  /(V'„)]/tan5„ 

and 

^n+i  =  On  -  2arctan/'(V'n)- 

Note  that  these  formulas  hold  if  x^  <  V’n  <  Xn+i-  Similar  formulas  can  be  derived  if  i„, 
V'n  and  in+i  are  permuted.  It  can  be  shown  that  (z„,  cos  6n)  are  canonically  conjugate  co¬ 
ordinates,  i.e.  the  one-point  ray  tracing  system  can  be  written  in  terms  of  a  Hamiltonian 
system. 

For  simplicity  we  took  a  simple  sinusoid  for  /,  with  a  wavelength  of  100  km  and  an  ampli¬ 
tude  of  0.1  km.  The  average  thickness  is  20  km,  the  same  as  used  in  the  layer  over  halfspace 
model.  Figure  5  shows  the  results  of  the  one-point  ray  tracing  for  SmS  through  7xSmS.  For 
SmS  the  first  multipathing  occurs  for  epicentral  distances  of  400  km. 

For  epicentral  distances  of  800  km  and  larger  it  becomes  hard  to  do  two-point  ray  tracing 
for  SmS.  If  the  undulation  of  the  Moho  becomes  larger,  then  it  turns  out  that  it  is  hard  to 
do  any  two-point  ray  tracing,  and  in  fact  the  ray  behavior  is  chaotic. 

Graph-theoretical  Rays  -  Because  of  these  difficulties  with  two-point  ray  tracing  in  media 
with  curved  interfaces,  it  is  worthwhile  to  investigate  the  possibility  to  use  other  methods, 
such  as  graph  theory.  Graph  theory  enables  fast  calculation  of  raypaths  and  travel  times  in 
comphcated  media,  avoiding  the  problems  mentioned  above,  and  makes  first  order  scattering 
computation  possible.  The  method  entails  approximating  ray  paths  with  shortest-time  paths 
between  nodes  in  a  grid  and  using  Dijkstra’s  algorithm  (Dijkstra,  1959)  to  search  over  all 
permissible  raypaths  between  a  source  and  a  receiver.  The  search  is  made  more  efficient, 
by  applying  a  sorting  algorithm  (Moser,  1991),  making  computation-time  proportional  to 
n  log  n,  where  n  is  the  number  of  nodes  in  the  model. 

The  method  finds  only  the  shortest-time  path;  ray  paths  representing  later  arrivals  can 
however  be  found  by  formulating  them  as  constrained  shortest  paths.  Scattered  ray-paths 


(2) 

(3) 
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can  be  found  by  requiring  the  rays  to  visit  a  set  of  nodes.  Moho  reflections  can  also  be  ob¬ 
tained  by  stacking  duplicates  of  the  crust,  the  number  of  duplicates  depending  on  the  number 
of  Moho  reflections  desired.  For  example,  to  obtain  ray  paths  and  travel  times  for  SmS  one 
duplicate  of  the  crust  is  required,  whereas  for  2x5m5  a  quadruplicated  stack  is  needed.  For 
one  reflection  the  model  consists  of  three  interfaces:  surface,  Moho  surface.  The  source 
is  located  at  a  node  on  the  surface,  on  one  side  of  the  model  and  the  shortest-time  ray  paths 
are  found  to  the  nodes  on  the  other  surface,  requiring  the  paths  to  go  through  the  nodes 
on  the  intermediate  interface,  representing  the  Moho.  We  have  made  such  an  adaptation  of 
the  method,  to  obtain  first  arrivals  for  multiple  Moho  reflections  of  SH  waves.  Our  initial 
model  is  the  same  ais  before.  The  nodes  are  1  km  apart  on  each  copy  of  the  Earth’s  surface 
and  the  Moho.  To  limit  the  ray  paths  to  Moho  reflections,  nodes  on  one  interface  are  only 
connected  to  nodes  on  the  next  one. 

Ray  paths  of  SmS  first-arrivals  in  this  model  are  shown  in  Figure  6.  As  distance  increases 
beyond  140  km,  the  bounce-points  of  the  first-arriving  rays  center  around  the  peaks  on  the 
Moho,  with  fewer  and  fewer  rays  leaving  the  source  ending  up  as  first  arrivals.  Thus  only 
a  fraction  of  the  energy  leaving  the  source  is  present  in  these  first  arrivals.  As  the  number 
of  Moho  multiples  increases,  even  less  of  the  source’s  energy  is  present  in  the  first  arrivals. 
Except  for  the  reflection  from  each  Moho  peak,  the  first  arrivals  in  this  model  are  not  geo¬ 
metric  rays,  as  the  reflection  angles  from  the  Moho  deviate  up  to  0.07  deg  from  Snell’s  law, 
but  rather  represent  diffracted  waves  from  the  peaks  on  the  undulating  Moho.  Travel-time 
deviations  from  those  of  a  flat  Moho,  are  extremely  small,  due  to  the  small  amplitude  of 
the  Moho  perturbation.  Therefore,  although  the  triplications  produced  by  the  sinusoidal 
Moho  are  not  obtained  with  this  method,  the  travel-time  curve  is  nearly  identical  to  the  one 
obtained  with  the  normal  ray  tracing  method. 

CONCLUSIONS  AND  RECOMMENDATIONS 

We  have  come  to  a  good  understanding  of  Lg-waves  from  both  a  mode  and  a  ray  point  of 
view  in  a  simple  layer  over  halfspace  model,  with  a  flat  discontinuity.  Prom  a  mode  point  of 
view  the  seismogram  consists  of  the  superposition  of  all  modes.  An  interpretation  in  terms  of 
Airy  phases  does  not  appear  to  be  useful.  From  a  ray  point  of  view  the  seismogram  consists 
of  all  crustal  multiples  that  are  critically  reflected.  The  agreement  between  the  seismograms 
obtained  using  these  two  methods  is  striking.  The  main  difference  is  that  head-waves  are 
not  modelled  by  the  ray  method  in  plane  layered  media. 

fr  the  Moho  is  undulating  then  the  two-point  ray  tracing  becomes  much  more  complicated. 
It  is  nevertheless  possible  to  do  two-point  ray  tracing  as  long  as  neither  the  amplitude  of  the 
undulation  nor  the  cpicentral  distance  arc  too  large.  The  graph  theoretical  traveltime  results 
are  nearly  identical  to  those  obtained  with  normal  ray  theory,  thus  providing  a  promising 
application  of  graph  theory. 

Further  research  should  focus  on  the  calculation  of  seismograms  in  media  with  one  or  pos¬ 
sibly  more  curved  interfaces.  The  possibilities  and  limitations  of  classical  ray  theory  should 
be  investigated,  as  well  as  the  supplementary  role  that  graph  theory  can  play.  It  is  also 
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necessary  to  compute  seismograms  in  more  complicated  layered  media,  using  both  methods. 
As  a  next  step  P-SV  waves  need  to  be  considered. 
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Figure  1.  Mode-seismograms  using  (from  top  to  bottom)  10 
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Figure  2.  Comparison  between  seismograms  computed  using  mode  theory  (solid)  and  ray  theory 
(dashed).  Epicentral  distances  are  100  km  (bottom),  300  km,  500  km  and  700  km  (top).  The 
reduced  velocity  is  the  S-wave  velocity  of  the  crust  (3.5  km/s). 
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.V  3  Comparison  between  seismograms  computed  using  mode  theory  (solid  lines)  and  ray 

theory  (daahed  lines).  The  epicentral  distance  is  500  km. 
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FigTire  4.  Geometry  of  layer  over  halfspace  model  with  undulating  discontinuity. 


Figure  5.  EpicentraJ  distance  as  a  function  of  take-off  angle  for  SmS  (bottom  curve)  through 
7x5m5  (top  curve).  Also  plotted  are  the  background  values  for  a  crust  with  a  thickness  of  20  km. 


Figure  6.  Ray  paths  for  SmS  from  a  surface  source  in  a  20-km-thick  crust,  with  sinusoidally 
undulating  Moho  of  amplitude  0.1  km  and  period  60  km.  Vertical  exaggeration  xlO. 
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ABSTRACT 

We  use  ray  theory  to  model  the  propagation  of  Lg  waves  through  2D  and  3D 
layered  crustal  models.  The  layers  are  homogeneous,  and  the  discontinuities  are 
undulating.  The  Lg  wave  train  is  modelled  by  multiple  S  reflections  within  the 
crustal  layers.  The  ray  tracing  system  is  reduced  from  a  set  of  linear  differential 
equations  to  a  set  of  maps.  If  the  medium  has  three  or  more  discontinuities 
the  number  of  multiples  increases  exponentially.  Using  a  binary  tree-searching 
method  we  systematically  keep  track  of  all  multiples. 

The  ray  behavior  is  chaotic  for  large  take-off  angles,  causing  multipathing 
of  the  crustal  multiples.  However,  in  the  presence  of  mid-crustal  discontinuities 
the  rms  Lg  amplitude  is  stable,  because  of  the  distribution  of  the  energy  over 
a  large  amount  of  rays.  The  simplest  case,  a  one  layer  over  a  halfspace  model, 
already  contains  all  the  characteristics  of  wave  propagation  through  the  media 
under  consideration. 

To  examine  the  source-depth  and  velocity-structure  effects  on  amplitudes 
and  character  of  Lg  we  use  wavenumber  integration  methods  in  1-D  structures. 
In  particular,  we  study  the  ability  of  shallow  explosions  to  generate  S*  and 
pS  rays.  Although  S*  is  a  nongeometrical  phase  it  does  follow  a  well  defined 
raypath.  In  velocity  structures,  where  the  peak  in  the  5*  radiation  pattern  is 
confined  to  the  crust,  signals  from  explosions  within  1  km  of  the  surface,  that  do 
not  generate  spall  or  significant  tectonic  release,  are  likely  to  be  dominated  by 
waves  with  a  very  narrow  ray-parameter  range,  coming  from  the  radiation  peak 
of  S*.  Observation  of  Lg  waves  dominated  by  such  a  narrow  ray-parameter 
range  may  therefore  identify  the  source  as  explosive,  and  shallow.  Dislocation 
radiation  patterns,  or  the  8p>allB  CLVD  term  are  likely  to  generate  significant  S 
waves  over  a  wider  range  of  ray  parameters  and  therefore  not  to  be  confused 
with  5*. 
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Effects  of  Explosion  Depth  and  Crustal  Heterogeneity  on  Lg  Waves 


OBJECTIVE 

The  objective  of  our  research  is  to  understand  aspects  of  the  excitation  of  Lg  waves  and  their  subsequent 
propagation  through  the  crust,  including  2D  and  3D  heterogeneous  crustal  models. 

Modeling  the  Lg  wave  train  is  a  difficult  task,  even  though  measurements  of  rms  amplitude  of  Lg  is  stable. 
We  attempt  to  explain  this  by  modeling  Xp  as  a  superposition  of  supercritically  reflected  crustal  multiples. 
In  our  models  the  crust  is  multilayered  (2D  or  3D)  and  the  interfaces  are  undulating. 

We  also  examine  the  effect  S*  has  on  amplitudes  and  character  of  Lg  to  determine  to  what  depths  and  in 
which  types  of  velocity  structures  S*  is  a  significant  contributor  to  Lg. 

RESEARCH  ACCOMPLISHED 

Lg  Wave  Propagation  Through  3D  Heterogeneous  Structures 

To  study  the  propagation  of  Lg  waves  through  2D  and  3D  heterogeneous  media  we  use  ray  theory.  Ray 
theory  is  not  only  very  efficient,  it  also  makes  the  characterisation  of  Lg  waves  in  terms  of  focusing/defocusing 
effects,  multipathing  and  traveltimes  due  to  the  complicated  structure  relatively  simple.  Thus  ray  theory 
complements  full  waveform  methods  such  as  finite  differences. 

Lg  waves  possess  a  dual  character  (e.g.  Hansen  et  al.,  1990).  On  one  hand  the  modeling  of  Lg  waves 
has  proved  to  be  very  difficult,  on  the  other  band  the  Lg  wavetrsdn  has  been  useful  to  determine  seismic 
moments  because  rms  Xj  is  a  stable  quantity,  that  is  relatively  independent  of  regional  structure.  Using  our 
modeling  results  we  give  an  explanation  for  this  seemingly  contradictory  character. 

The  models  that  we  use  are  2D  and  3D  layered  crustal  models.  The  layers  are  homogeneous,  and  the 
discontjnujtiw  are  undulating.  The  Lg  wave  train  u  modelled  by  multiple  5  reflections  within  the  crustal 
layers.  In  a  flat  layered  medium  the  ray-theoretical  seismograms  are  almost  identical  to  the  exact  synthetics 
(Keers  et  al..  1995a).  This  serves  as  our  justification  for  the  use  of  ray  theory  in  media  with  (slightly) 
UDdul&tiDg  iDtcrf&ccB  •cp&iAtcd  by  homogCDcouB  Iftycn. 

One  of  the  mam  advantages  of  using  homogeneous  layers  is  that  the  ray  tracing  system  can  be  reduced 
from  a  set  of  linear  differential  equations  to  a  set  of  maps  (Keers  et  al.,  1994).  This  makes  the  ray  tracing 
very  efficient;  it  is  not  necessary  to  use  a  Runge-KutU  method  and  the  search  for  the  point  where  the  ray 
intersects  one  of  the  undulating  interfaces  u  reduced  to  a  ID  root-solving  problem  (even  in  the  three  dimen¬ 
sional  case).  If  the  medium  has  three  or  more  discontinuities  the  number  of  multiples  increases  exponentially 
Using  a  binary  tree-searching  method  we  systematically  keep  track  of  all  multiples. 

The  simplest  case,  a  one  layer  over  a  halfspace  model,  already  contains  all  the  characteristics  of  wave  prop¬ 
agation  through  the  media  under  consideration.  In  general  two  regions  (by  regions  we  mean  take-off  angle 
intervals)  are  distinguished.  One  region,  corresponding  to  smaller  take-off  angles,  has  regular  ray-behavior 
and  another  region,  corresponding  to  the  largest  take-off  angles,  has  chaotic  ray-behavior.  In  the  case  of 
a  sinusoidal  Moho  with  a  wavelength  of  100  km  and  an  amplitude  of  1  km  the  take-off  angle  intervals  arc 
roughly  50»-65»  and  65».90»  (the  critical  angle  U  close  to  50»).  The  first  region  is  characterized  by  a  low 
degree  of  multipathing  and  a  high  degree  of  focusing.  The  second  region  is  characterized  by  a  large  degree  of 
rnultipathmg  and  consequently  a  large  degree  of  defocusing  The  defocusing  is  caused  by  the  fact  that  a  small 
change  in  take-off  angle  causes  a  large  change  in  epicentral  distance,  except  at  the  caustics.  However,  at  the 
caustics  the  amplitude  is  still  relatively  small  (Keers  et  al,  1995a).  It  should  be  noted  that  the  amplitude 
perturbations  due  to  the  undulation  are  much  larger  than  the  traveltime  perturbations.  Table  1  gives  an 
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indication  of  the  degree  of  multipathing  for  the  case  of  a  flat  Moho,  a  sinusoidal  Moho  and  a  more  realistic 
model,  that  of  the  Moho  below  Germany. 

In  the  case  of  a  more  realistic  two-layered  model  with  undulating  interfaces,  the  ray  behavior  remains  es¬ 
sentially  the  same:  there  exists  a  regular  region  corresponding  to  smaller  take-off  angles  and  a  chaotic  region 
corresponding  to  larger  take-off  angles.  All  discontinuities  are  sinusoidal.  The  surface  has  a  wavelength  of  60 
km  and  amplitude  of  0.5  km,  the  midcrustal  discontinuity  has  a  wavelength  of  160  km  and  an  amplitude  of 
0.8  km  and  the  Moho  has  a  wavelength  of  100  km  and  an  amplitude  of  1.5  km.  Figure  la  shows  the  amplitude 
distribution.  This  figure  is  essentially  an  ’envelope  record-section’.  It  contains  only  amplitude  information; 
no  phase  information.  The  strong  focusing  in  the  regular  region  that  was  present  in  the  one-layered  model 
is  gone.  This  is  due  to  the  reflection  and  transmission  of  the  multiples  at  the  midcrustal  interfaces,  which 
partitions  the  energy  over  many  multiples.  The  presence  of  the  chaotic  region  makes  it  impossible  to  identify 
even  the  major  multiples  (SmS,  2SmS  etc.)  at  larger  epicentral  distances.  The  large  degree  of  multipathing 
of  a  certain  multiple,  like  2SmS,  causes  these  multiples  to  consist  of  many  rays,  ejirli  with  a  relatively  small 
amplitude,  coming  in  from  totally  different  directions.  Each  of  these  multiples  is  smeared  out  over  a  larger 
time  interval  and  the  the  maximum  amplitude  of  the  dominant  rays  is  low  compMued  to  the  flat  layered 
case.  This  explains  why  the  modeling  of  Lg  waveforms  has  been  largely  unsuccesful:  the  waveform  is  very 
sensitive  to  initial  conditions  (take-off  angle)  and  model  parameters. 

Figure  2a  shows  the  effect  of  a  3D  sinusoidal  Moho.  The  wavelength  is  100  km  in  both  the  x  and  y 
directions,  and  the  amplitude  is  1  km.  The  figure  shows  the  energy  distribution  along  the  x  axis.  Also 
shown  is  the  energy  distribution  for  a  2D  model  with  a  sinusoidal  Moho  that  has  a  wavelength  of  100  km 
and  amplitude  of  1  km.  The  amount  of  side  scattering  can  not  be  ignored,  but  depends  strongly  on  the 
model  (Keers  et  al,  1095b).  Since  this  3D  example  is  for  a  one-layered  model,  the  focusing  of  the  different 
multiples  plays  a  dominant  role  in  this  figure. 

Figure  2b  shows  the  total  energy  of  the  Lg  wavetrain  as  a  function  of  epicentral  distance  for  the  two 
models  of  figure  1.  Also  plotted  is  the  energy  of  the  waves  assuming  a  (Nuttli,  1973)  decay  of  the 

amplitude  (A=distance).  The  undulating  twc^layered  model  shows  much  less  focusing  especially  at  larger 
epicentral  distances.  This  is  due  to  the  midcrustal  interface.  There  is  still  focusing  of  certain  rays,  but  the 
reflection  and  transmission  from  the  midcrustal  interface  makes  these  rays  much  less  dominant  than  in  the 
one-layered  model.  This  explains  the  sUbility  of  rms  Lg  as  observed  for  many  regional  events  (e.g.  Hansen 
et  al.,  1990).  For  models  with  more  layers  we  can  expect  the  focusing  effects  to  become  even  less  dominant. 

Effects  of  Source  Depth  and  Velocity  Structure  on  tbe  Character  of  Lg 

Discrimination  between  shallow  sources  is  a  difficult  task  and  methods,  such  as  Pg/Lg  amplitude  ratios 
and  Lg  spectral  ratios  as  well  as  regional  seem  to  be  region  dependent  and  only  work  for  discrimina¬ 

tion  between  shallow  explosions  and  earthquakes  deeper  than  5  km,  and  even  then  there  are  problems  with 
explosion  depths  of  <1  km.  Generation  of  S  waves  from  explosion  sources  has  generally  been  attributed  to 
nonisotropic  source  radiation,  spall,  tectonic  release,  or  fig-S  scattering  in  the  source  region,  while  S' ,  the 
wave  generated  by  interaction  between  the  curved  P  wave-front  and  the  surface,  has  usually  been  discounted 
for  all  but  the  shallowest  sourcs,  because  of  its  strong  dependence  on  source  depth.  However,  due  to  the  5* 
radiation  pattern,  in  certain  structures  the  amplitude  may  be  significant  down  to  depths  of  ~1  km,  possibly 
making  S*  an  added  tool  to  discriminate  between  shallow  sources. 

To  examine  the  source-depth  and  velocity-structure  effects  on  amplitudes  and  character  of  Lg  we  use 
wavenumber  integration  methods  in  1-D  structures.  The  shortcomings  of  the  simplified  structure  are  over¬ 
come  by  the  benefiu  of  a  full-wave  synthesis,  revealing  all  waves  generated  by  the  source  and  its  interaction 
with  discontinuities  in  the  structure;  of  main  interest,  of  course,  u  the  interaction  with  the  Earth’s  surface. 
The  characteristics  observed  in  the  synthetics  are  then  explained  in  terms  of  the  radiation  pattern  of  S’. 

A  high  frequency  approximation  to  the  amplitude  and  phase  of  the  S’  wave  was  developed  by  Daley  and 
Hron  (1983).  In  this  approximation  5*  appears  as  a  spherical  wave-front  with  5V-particle  motion,  radiated 
from  a  point  on  the  surface  above  the  source.  lU  amplitude  is  modulated  as  a  function  of  angle  d  (measured 
from  the  vertical),  but  is  independent  of  aiimuth.  The  wave  exists  for  4>i  =  sin“^(/3/a)  <  4>  <  r/2,  but 
the  approximation  is  valid  for  (t>i  <  <!><  t/2,  or  l/o  <  p  <  1//3,  where  a  and  ^  are  the  m^iums  P  and  S 
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velocities  and  p  =  sin  is  the  ray  parameter.  A  plot  of  the  radiation  function  for  a  frequency  of  2  H* 
and  a  =  4.2  km/s  is  shown  in  Figure  3  (top  right).  The  source  depths  represented  are:  0.1,  0.5,  1.0  and  2.0 
km.  The  amplitude  is  a  maximum  at  but  decreases  to  sero  at  ^  =  45®,  where  the  phase  also  flips  by  90®. 

The  amplitude  increases  agun  to  a  smaller  maxi  mum,  before  becoming  sero  again  at  90®.  The  amplitude 
exponentially  with  source  depht,  so  for  depths  greater  than  0.5  km,  amplitudes  are  very  small  at 
angles  above  45®.  However,  between  and  45®  the  amplitude  can  be  larger  than  in  the  original  P  and  the 
reflected  pS  wave.  To  demonstrate,  the  pS  radiation  pattern,  which  exists  for  reflection  angles  ^  is 

also  included  on  the  plot. 

If  the  narrow  radiation  peak  of  S*  gets  trapped  in  the  crust,  it  may  dominate  Lg.  To  demonstrate  we 
calculated  t  juthetics  in  three  diflerent  velocity  structures,  mdl,  where  the  peak  is  radiated  irito  the  mantle, 
md2,  where  the  peak  is  confined  to  SmS  multiples  for  source  depths  down  to  1.25  km  and,  md3,  where  the 
peak  is  confined  to  turning  waves  in  the  crust  for  source  depths  down  to  1  km.  The  models  have  a  common 
Q  structure  and  share  a  velocity  structure  below  3  km  depth  (see  Figure  4,  lower  panel).  It  is  the  ratio 
Q*erwre«/Anonfie  that  determines  how  much  of  5*  is  trapped  in  the  crust.  If  a««urcc  <  Pmantu>  ^  of  S'  and 
some  pS  arc  contained  in  Lg,  but  when  Ojourc*  increases  to  >  Pmantu  the  radiation  peak  starts  to  go  into 
the  mantle.  When  the  critical  angle  for  Moho  reflection  exceedes  ~  42®,  no  significant  S'  can  be  found  in 
Lg,  except  for  depths  <0.6  km.  The  angular  range  of  Moho  reflections  in  the  three  models  is  indicated  on 
Figure  3  (top  right).  There  it  is  clear  that  in  mdl,  most  of  the  S'  radiation  peak  goes  into  the  mantle  smd 
therefore  Lg,  from  a  pure  explosion  source,  is  expected  to  decay  quickly  with  source  depth.  However  in  md2, 
all  of  the  radiation  peak  ends  up  in  Moho  reflections,  for  source  depths  down  to  1.25  km.  For  this  depth 
range  Lg  is  expected  to  be  of  significant  amplitude  and  to  be  dominated  by  Moho  reflections.  In  md3,  the 
peak  is  confined  to  turning  waves  in  the  crust,  for  source  depths  down  to  1  km.  In  this  depth  range  Lg  is  also 
expected  to  be  significant  and  dominated  by  turning  waves.  These  effects  are  clearly  seen  in  the  synthetics: 
In  Figure  4,  synthetics  are  shown  at  a  distance  of  500  km  for  source-depths  ranging  from  0.1  to  2  km  in 
model  md2.  Cle^ly  the  Lg  wave  is  dominated  by  the  Moho  multiples,  2SmS,  ZSmS  and  4SmS,  and  Lg 
amplitude  is  significant  for  source  depths  down  to  ~1  km.  Record  sections  for  an  explosion  source  at  0.5  km 
depth  in  the  three  models  are  shown  in  Figure  5,  where  the  expected  features  are  also  seen.  The  mdl  record 
section  has  smaller  amplitude  Lg  waves  than  the  other  models,  evenly  divided  between  Moho  reflections  and 
turning  waves.  The  frequency  content  of  the  Lg  wave  it  also  significantly  lower  than  that  of  the  crustal  P 
wave-train.  The  md2  record  section  has  much  larger  amplitudes,  Lg  is  dominated  by  Moho  multiples,  and 
there  is  no  significant  loss  of  higher  frequencies  in  Lg.  The  md3  record  section  has  large  amplitude  Lg  waves 
of  high  frequency  content,  and  maximum  amplitudes  coincide  with  the  arrivals  of  the  turning  waves  and 
Moho  reflections;  the  Moho  reflections  come  from  pS  and  the  turning  waves  come  from  S'. 

Analysis  of  Array  Data 

The  theoretical  results  need  to  be  tested  against  observed  data  from  local  and  regional  events.  We  have 
so  far  selected  two  regions  based  on  the  availability  of  short-period  array  data.  The  regions  are:  Central 
Europe,  where  we  make  use  of  the  GERESS  array,  and  the  Sierra  Nevada,  where  three,  3-component  arrays 
were  temporarily  operated  in  1993.  The  arrays  enable  detection  of  small  coherent  arrivals  and  decomposition 
of  the  Lg  wave  train  into  distinct  arrivals  with  differing  phase  velocities.  The  events  recorded  at  GERESS 
consist  of  earthquakes  and  explosions,  while  most  of  the  events  recorded  on  the  Sierra  Nevada  arrays  are 
earthquakes,  but  also  include  the  NPE  explosion  at  the  Nevada  Test  Site  in  September  1993. 

CONCLUSIONS  AND  RECOMMENDATIONS 

In  2D  and  3D  heterogeneous  media  Lg  wave  propagation  is  strongly  affected  by  curved  interfaces.  The 
ray  behavior  is  chaotic  for  large  take-off  angles,  causing  multipathing  of  the  crustal  multiples.  This  makes 
modeling  of  Lg  waves  a  tidy  task.  However,  in  the  presence  of  mid-crustal  discontinuities  rms  Lg  is  stable, 
because  of  the  distribution  of  the  energy  over  a  large  amount  of  rays.  It  should  be  noted  that  the  undulations 
used  are  smaU  (of  the  order  of  1  km).  The  amplitude  behavior  wUl  change  if  the  undulations  of  the  Moho 
become  large.  Further  research  is  needed  to  establish  this  effect.  Incorporation  of  P-S  and  S-P  conversions 
using  our  ray-tracing  algorithm  should  make  it  possible  to  model  Lg  wave  propagation  more  realistically. 
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More  de^ed  Moho  mapt  for  regioni  of  ioterot  are  al^  neceaMry  for  reliable  modeling  of  Lg  in  3D  atruc- 
tum.  This  goal  can  be  achieved  by  wave-form  inversion,  using  regional  seiamograms  (Das  and  Nolet,  1995) 
In  velocity  rtrurt,^.  where  the  peak  in  the  5*  radiation  pattern  is  confined  to  the  cniat,  aignaJs  from 
^^lons  wit^  1  km  of  the  «irface,  that  do  not  generate  apaU  or  significant  tectonic  relea«:,^  likely 
to^  donated  by  wava  with  a  very  narrow  ray-parameter  range,  coining  from  the  radiation  peak  of  5* 
Observation  of  Lg  wara  donated  by  such  a  narrow  ray-parameter  range  may  therefore  identify  the  source 
as  «^lo»ive.  Md  shallow.  Dislocation  radiation  patterns,  or  the  spalls  CLVD  term  are  likely  to  generate 
significant  S  waves  over  a  wider  range  of  ray  parameters  and  therefore  not  to  be  confused  with  S*  A 
comparison  between  daU  from  mining  explosions  and  earthquakes,  located  m  approximately  the  same  area 

may  serve  as  a  c«e  for  this.  GERESS  daU  for  events  in  the  Vogtland  region  in  the  Csech  Republic 
may  be  powible  CAndidates. 

The  completed  record  section  profiles  will  serve  as  the  data  base  needed  to  test  the  validity  of  the  theo¬ 
retical  calculations  and  predictions. 
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Table  1 

Number  of  supercritically  reflected  rays  for 
different  Moho  models 
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Figure  3.  (upper  right)  5*  radiation  pattern  at  2  Hr,  for  varioui  iourcc  dcpthi  in  a  medium  with  o  =  4.2  km/t. 
Velocity  modeit  (left)  and  Q  structure  (lower  right)  used  in  the  synthetic  calculations.  Stars,  representing 
some  of  the  source  dcpthi  used,  are  positioned  at  the  sub-Moho  5  velocity,  Pmmntu- 


Figure  4.  Vertical  component  lynlhetic  duplacement  at  500  km  distance  from  an  explosion  source  at  various 
source  depths  in  model  md2.  The  synthetics  have  a  Nyquist  frequency  of  4  Hs  and  are  convolved  with  a 
NORESS-type  shorUperiod  instrument. 


ffldl 


fl  -  M/7A]  imm) 


r  -  i/TJ  )  imm) 
m4^ 


Figure  5.  Venice]  compoDcnt  lynthctic  record  »cction«  for  an  explosion  depth  of  0.5  km  in  the  three  models, 
mdl,  md2  ajid  md3,  with  tr&vcl-timc  curvet  for  the  major  phases  superimposed.  Amplitudes  are  normalized 
to  their  maxjma,  shown  al  the  right  of  each  trace.  Slopes  on  the  travel-time  curves  can  be  read  from  the 
velocity  template  in  the  lower  right  corner. 
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ABSTRACT 

This  report  discusses  both  theoretical  and  observational  aspects  of  asymptotic  crustal  wave  prop¬ 
agation. 

Crustal  models  (both  2D  and  3D)  are  conveniently  expressed  using  undulating  interfaces.  This 
makes  it  necessary  to  develop  special  ray  tracing  algorithms  (both  for  the  kinematic  and  dynamic 
ray  tracing).  These  algorithms  and  the  chaotic  properties  of  rays  traveling  through  these  models 
are  discussed  briefly. 

One  consequence  of  the  chaotic  ray  behavior  is  the  large  number  of  caustics  that  are  generated. 
This  makes  it  necessary  to  use  Maslov  theory  to  compute  synthetics,  rather  than  geometrical  ray 
theory.  It  is  shown  in  this  report  how  to  compute  the  correct  Maslov  phase  and  amplitude  using 
the  method  of  stationary  phase. 

One  of  the  major  obstacles  in  computing  Maslov  synthetics  has  been  the  presence  of  spurious 
arrivals  caused  by  pseudo-caustics.  Using  the  method  of  stationary  phase  these  pseudo-caustics 
are  shown  to  disappear  if  one  uses  2D  Maslov  integrals.  An  algorithm  is  presented  that  efficiently 
computes  2D  Maslov  integrals.  Previously  Maslov  synthetics  were  plagued  by  the  effect  of  the 
endpoints  of  the  integration.  This  endpoint  effect  is  removed  by  a  simple  change  of  variables. 

The  Maslov  synthetics  can  be  used  to  justify  the  use  of  ray  density  plots  as  energy  plots.  The 
ray  density  plots  are  used  to  study  array  data  from  the  NPE  at  the  Nevada  test  site.  Lg  waves, 
propagating  in  the  direction  of  the  southern  Sierra  Nevada  range,  disappear  in  the  epicentral 
distance  r  mge  of  250  km  -  350  km.  It  is  shown  how  the  theoretical  modeling  of  Lg  can  be  used  to 
put  constraints  on  the  crustal  structure  of  sourthem  California. 

Keywords:  crustal  waves,  chaotic  rays,  Maslov  synthetics,  energy  plots,  array  data 
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OBJECTIVES 


The  first  goal  of  this  project  is  to  develop  reliable  and  efiicient  methods  to  compute  synthetic 
seismograms  through  2D  and  3D  heterogeneous  crustal  models,  using  asymptotic  methods.  The 
second  goal  of  this  project  is  to  apply  these  methods  to  study  crustal  wave  propagation  (with 
special  emphasis  on  Lg)  using  array  data.  Both  of  these  aspects  are  discussed  in  this  report. 

RESEARCH  Av.€OMPLISHEI> 


Kinematic  Ray  Tracing  Through  3D  Multilayered  Media 

We  briefly  describe  raytracing  in  3D  layered  media  with  undulating  interfaces  separated  by  homo¬ 
geneous  layers.  Two  aspects  are  important;  first  one  needs  to  trace  rays  from  one  interface  to  the 
next  interface.  Secondly  all  multiples  generated  need  to  be  traced.  Here,  for  the  sake  of  conciseness, 
we  only  consider  the  first  aspect  of  the  algorithm;  for  more  details  see  Keers  et  al.  (1996b). 

The  interfaces  are  described  by  the  equation  z  =  hp  +  fp{x,y)  and,  as  before,  the  pth  layer  is 
right  below  the  pth  interface.  We  choose  the  coordinate  system  as  shown  in  figure  1. 

Consider  a  ray  impinging  upon  the  pth  interface  from  above.  The  unit  vector  pointing  in  the 
direction  of  the  rayisin  =  (cos  0:^1 » cos cos  7„)  (see  figure  1).  This  vector  is  assumed  to  be  known. 
The  unit  normal  to  the  pth  interface  at  (i„, p„,  /ip  +  /p(i„,  p„))  is  denoted  by  n  =  (ni,  712,  ng).  Prom 
figure  2  it  can  be  seen  that  r„  =  2(i„  •  n)n  —  i„,  where  r„  is  the  unit  vector  pointing  in  the  direction 
of  the  reflected  ray.  The  derivation  of  t„  is  a  bit  more  involved.  Prom  figure  2  we  see  that 


t„  -  (t„  •  n)n  =  c(i„  -  (i„  •  n)n),  (1) 

for  some  c  <  0.  We  compute  c  by  taking  the  inner  product  of  (1)  with  in  and  t„  respectively.  We 
find 

tn  =  (tn  •  n)n  +  t„  —  (t„  •  n)n  =  —  ^  ^  ^  }  (^n  '  ~  cos^^  n.  (2) 

Here  “tp  is  found  with  Snell  s  law.  These  results  enable  us  to  derive  the  raytracing  equations.  If  a 
ray  goes  up  then  direction  vector  is  denoted  by  u^.  If  the  ray  goes  down  then  the  direction  vector 
is  d^.  It  can  be  shown  that 


^  fpi^niVn)  + 


(!:)<» 


p-1 

n+1 


vF  \ 


■n+l 


O)  (3) 


which  is  an  implicit  equation  for  ,and 


-  =  y^rTA- 


(4) 


The  reflection  vector  d^.,.1  and  the  transmission  vector  are  found  using  the  formulas  given 
above  and  the  correct  choice  of  orientation  of  n.  In  vector  notation: 


(6) 


dn+1  =  2(uP  •  n)n  +  < 
ap-2 


uK  =  (^)  +  (  (^)  «  •  n)  -  |1  -  (^)  \l  -  «  .  n)^)lV^ j  n. 

Similar  formulas  can  be  derived  for  the  downgoing  ray  (Keers  et  al,  1996b). 
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Dynamic  Ray  Tracing 

We  now  proceed  to  derive  the  partial  derivatives  of  the  kinematic  ray  equations.  It  is  clear  that  in 
3D  media  the  parameters  that  give  the  ray  positions  are:  i„,  and  d„  (the  superscripts  p  and 
p  -  1  are  deleted  from  now  on).  AH  the  other  parcimeters  are  found  from  these.  In  particular  we 
shall  be  interested  in  the  dependence  of  these  parameters  on  the  initial  take-off  angle  Bq.  These  are 
necessary  in  establishing  the  degree  of  chaoticity  using  Lyapunov  exponents  and  the  computation  of 
synthetics  with  Maslov  integrals  (see  the  next  two  sections).  In  3D  media  dn  has  three  components 
and  for  the  computation  of  the  Maslov  amplitude  we  also  need  the  dependence  of  i„,  and  d„  on 
the  initial  azimuth  ^o-  The  differentation  can  be  computed  by  perturbing  the  ray  equations  with 
the  appropriate  parameter  (e.g.  8xn)  and  then  derive  expressions  of  the  form 


^^n  +  l 


Taking  the  limit  gives  the  required  partial  derivative.  For  brevity  we  only  give  some  of  the  partial 
derivatives.  Note  that,  whenever  possible,  use  should  be  made  of  the  vector  notation.  As  in  the 
case  of  the  kinematic  ray  tracing  this  considerably  reduces  the  algebraic  effort.  One  finds 


_  djp-l  gVn^l  ^  ^  ,  a/p-l  dpn+l  \ 

dxn  V  ^  dy„+i  dxn  J  dyn+i  dx^  )  ’ 
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etc.  where 


dXn 


2X“t,  5i,n)n  -r  (u„  ■  n)di„n]. 


etc.  Similar  partied  derivatives  are  derived  for  the  upgoing  ray. 


(10) 

(11) 

(12) 


Chaotic  Ray  Behavior  and  Lyapunov  Exponents 

We  showed  previously  (Keers  et  al.,  1996a)  that  a  small  change  in  the  take  off  angle  gives  rise  to  a 
large  change  in  the  epicentral  distance.  This  chaotic  behavior  can  be  quantified  using  the  concept 
of  Lyapunov  exponents.  The  idea  is  to  measure  the  separation  d  between  two  points  in  phase  space 
that  are  initially  together.  Note  that  two  points  in  phase  space  axe  close  together  only  if  they  are 
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at  almost  the  same  position  in  coordinate  space  and  if,  in  addition  to  that,  the  slownesses  of  the 
rays  emanating  from  these  points  are  almost  the  same  as  well. 

If  the  ray  tracing  system  is  chaotic  then,  by  definition,  the  separation  rate  of  the  two  points  in 
phase  space  is  exponential.  In  other  words  there  is  a  constant  I  >  0  such  that 


e"^  =  d(n),  (13) 

In  order  to  illustrate  this  we  use  the  2D  version  for  a  one  layered  model  of  the  kinematic  and 
dynamic  ray  equations  derived  above.  After  scaling  i  with  the  thickness  of  the  crust  one  can  define 

■^(^ni^n)  —  (^n+1  >  ^n+l)  B-nd 
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Let  —  {6n,Xn)  and  =  ((i0„,  S6„  and  ^iq  are  small  then  yo  and  yo  +  Syo  are  close 

together.  After  one  iteration  we  have  y^  =  Af  (yo)-  Since  |  Syo  |  is  very  small  M{yo  +  6yo)  may  be 
approximated  by  M{yo)  +  D{yo)6yo  =  yi  +  D{yo)6yo.  Therefore  6yi  a  D{yo)Syo.  More  generally 
^J/n  =  D(y„^i)6yn-i  =  D”'(yo)6yo,wh.eTe  Z)"(yo)  =  ^iyn-i)D{yn-2)'D{yo)-  The  distance  between 
yo  and  yo  +  6yo  after  n  iterations  can  therefore  be  approximated  by  |  ^y„  |=|  D"(yo)(5yo  |.  The  rate 
of  divergence  I  between  yo  and  yo  and  yo  +  6yo  is  therefore  given  by 


^ni  _  Uyn  I 

Uyol 


(15) 


In  principle  I  depends  on  yo,  6yo  and  n:  I  =  /(yo,  Syo,  n).  However  /  can  be  shown  to  be  independent 
of  the  initial  conditions,  and  to  converge  to  a  constant  as  n  -*  oo.  Therefore  it  is  useful  to  define 
the  Lyapunov  exponent  /£,: 


Il{0o,xo) 


1  L>"(yo)<?yo 
n  I  (5yo  i 


(16) 


When  explicitly  computing  the  Lyapunov  exponents  we  use  the  method  of  Benettin  et  al.  (1980). 
Choose  (SBofSxo)  arbitrarj'  and  rescale  {66,y,6xn)  every  p  iterations  by  dividing  by  its  magnitude. 
In  this  way  numerical  overflow  is  avoided.  If  we  do  this  q  times  and  denote  the  resulting  scaling 
factors  by  s,(t  =  then 


In  s,  . 


»=i 


(17) 


In  practice  we  choose  p  =  10  and  q  =  100.  Figure  3b  shows  the  results  of  the  computation  of  9 
Lyapunov  exponents.  The  corresponding  solutions  of  the  rays  in  phase  space  are  shown  in  figure 
3b.  The  Lyapunov  exponents  of  all  the  rays  converge  to  a  certain  value.  In  fact,  since  they  converge 
relatively  rapidly,  we  could  have  taken  q  =  50.  The  most  chaotic  ray  in  figure  3b  with  an  initial 
take-off  angle  of  70“ ,  has  a  Lyapunov  exponent  of  0.8.  This  means  that  rays  that  are  initially 
100  m  apart,  after  ten  reflections  from  the  Moho  are  almost  300  km  apart;  this  strongly  chaotic 
behavior  results  in  a  large  degree  of  multipathing  and  many  caustics.  When  computing  synthetics 
It  therefore  is  necessary  to  use  an  algorithm  that  easily  can  handle  these  phenomena.  The  major 
candidate  for  this  is  Maslov  theory. 
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Maslov  Theory  and  Pseudo-caustics 

An  acoustic  wavefield  u  excited  by  a  point  source  at  xq  at  time  0  satisfies 

V^u(x,t)-  =  S(x  -  xo)6(t)  (18) 

If  we  define  the  Fourier  transform  as 

/(u,)  =  r  f{t)e-^dt  (19) 

J  —  OO 

then  the  wave  equation  can  be  transformed  to  the  frequency  domain  and  we  have 

V^u(x,a;)  +  k^{x)u{yc,uj)  =  ^(x  -  Xq) 

Here  k  =  uj/c.  hi  Maslov  theory  one  looks  for  solutions  of  the  form 

u{x,u;)=  J  A{x,q)e'^'^^^'’^^dq  (21) 

and 

V.{x,U})  = 

Whenever  both  asymptotic  ray  theory  (ART)  and  one  of  the  Maislov  integrals  are  valid  then  the 
ART  solution  can  be  obtained  from  the  Maslov  integrals  by  applying  the  method  of  stationary 
phase.  One  advantage  of  this  is,  that  the  Maslov  amplitude  can  easily  be  found.  Application  of 
the  method  of  stationary  phase  to  the  ID  Maslov  integral  gives  the  expression  for  A.  If  we  take  q 
to  be  the  initial  take-off  angle  </>o,  then  the  wavefield  becomes 

u(x,t4^)  =  ^  A€^'^d(f>Q  (23) 

This  has  the  advantage  that  the  integration  is  over  a  closed  circle.  The  usual  bounda  ry  points 
therefore  don’t  exist,  and  there  are  no  spurious  arri^'als  in  the  integral.  Also  the  application  of  a 
taper  at  the  boundaries  becomes  superfluous.  The  stationary  phase  method  gives  A  again.  Also  in 
this  case  the  we  may  take  (91,^2)  =  (<^i^o)t  in  which  case  the  integral  reduces  to 

u{x,u;)  =  ^A{x,<Po,eo)e^'^^^'‘^-^°'>d4>odeo  (24) 

It  is  worthwile  to  have  a  closer  look  at  the  singularities  of  the  asymptotic  wavefield  expressions. 

Classical  ray  theory  breaks  down  at  points  where  J  =  0  ( J  is  the  Jacobian  from  coordinate  space 
to  ray  coordinates).  Similarly  the  Maslov  integrands  are  singular  at  points  where  J’  =  0  (here  S  is 
the  Jacobian  from  mixed  phase  space  to  ray  coordinates).  These  points  form  pseudo-caustics.  In 
the  case  of  a  ID  integral  these  points  also  correspond  to  stationary  points.  In  the  2D  integral  the 
pseudo-caustics  are  given  by 

\J\  =  0  (25) 


j  j  A{x,qi,q2)e*^^^^’'^'-’^'>dqidq2 
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However,  stationary  points  are  given  by, 

=  (26) 

||(.-^)+||(j-r)  =  0.  (27) 

A  singul'  rity  in  the  Maslov  integral  then  would  only  be  prominent  whenever  v?’e  have  simultaneously 
a  pseudo-caustic  and  a  stationary  point.  Therefore  pseudo-caustics  are  less  troublesome  in  2D 
integral  expressions  than  in  ID  integral  expressions. 

In  the  time  domain  we  get  the  following  expression  for  the  Maslov  wavefield 

u(x,  i)  =  ^Re[A(<) .  ^  A  ]  (28) 

where  A(t)  =  S(t)  -  i/rt  is  the  analytic  delta  function.  In  the  ID  case  the  integration  is  done 
similarly  to  the  WKBJ  evaluation  of  wavefields  in  ID  media  following  Chapman  et  al.  (1988). 
In  the  2D  case  we  either  are  integrating  over  a  plane  or  a  sphere,  as  discussed  previously.  This 
requires  a  slight  adaptation  of  the  integration  algorithm  (subroutine  THETAC  (Chapman  et  al., 
1988).  The  integration  domain  is  triangulated  (Delaunay  triangulation  for  example)  and  a  suitable 
taper  is  chosen  in  the  case  that  the  integration  domain  is  plane.  The  Delaunay  triangulation  has  the 
useful  property  that  the  3  gridpoints  closest  to  each  point  in  a  triangle  form  the  triangle.  The  best 
linear  approximation  to  a  function  at  a  point  in  the  triangle  that  uses  the  gridpoints,  is  therefore 
given  by  the  average  of  the  values  of  that  function  on  the  three  vertices  of  the  triangle. 

Ray  Density  Plots  and  Lg  Propagation  in  California 

The  Maslov  method  can  also  be  used  to  justify  the  use  of  ray  density  plots  as  energy  plots.  Energy 
plots  are  useful  for  studying  the  crustal  wavetrain  as  a  whole.  Figure  4  shows  vertical  component 
data  from  the  Non  Proliferation  Experiment  (NPE)  at  the  Nevada  Test  Site.  All  data  are  highpassed 
at  2  s.  Traveltime  curves  of  PmP  and  SmS  are  shown  as  well  as  traveltime  curves  of  other  dominant 
crustal  multiples.  The  source  and  station  locations  are  shown  in  figure  5.  We  concentrate  in  this 
report  on  Lg.  Using  a  three  layered  crustal  model  ray  density  plots  are  obtained  for  different  models. 
In  the  case  the  model  is  flat  the  different  multiples  can  be  readily  identified  (figure  6a).  However 
from  the  data  it  can  be  seen  that  Lg,  while  stiU  strong  at  270  km,  quickly  reduces  in  amplitude  and 
is  almost  mvisible  in  several  records  at  distances  larger  than  300  km.  More  realistic  models  include 
crustal  thickening  under  the  Sierra  Nevada  and  bumps  in  the  midcrustal  interfaces.  The  resulting 
ray  density  plot  for  a  10  km  amplitude  and  150  km  wavelength  bump  is  shown  in  figure  6b.  If,  in 
addition  to  these  undulations,  variations  in  the  midcrustal  interfaces  are  added  (figure  6c)  then'ug 
disappears  at  the  required  epicentral  distance.  Using  the  ray  density  plots  it  is  possible  to  constrain 
the  heterogeneity  of  the  Moho  between  the  NTS  and  the  receivers.  It  was  found  that  the  bump  in 
the  Moho  has  its  maximum  just  slightly  to  the  east  of  the  center  of  the  Sierra  Nevada.  This  result 
is  rather  independent  of  the  strength  of  the  heterogeneities  in  the  midcrust.  Figure  6d  shows  the 
ray  density  plots  for  a  model  like  that  in  figure  6c,  but  with  additional  random  undulations  of  the 
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midcrustal  interface  with  an  amplitude  of  2  km  and  and  a  dominant  wavelength  of  100  km.  The 
random  structure  has  the  effect  to  even  out  some  prominent  arrivals  still  present  in  figure  6c.  Since 
the  data  do  not  show  such  arrivals,  we  prefer  the  last  model. 

CONCLUSIONS  AND  RECOMMENDATIONS 

We  show  in  this  report  how  to  compute  ray  paths  through  2D  and  3D  layered  crustal  structures. 
These  crustal  models  generate  cha^otic  rays,  the  effect  of  which  is  to  create  many  caustics  and  a 
high  degree  of  multipathing.  Using  Maslov  theory  synthetics  can  be  computed;  the  effect  of  pseudo¬ 
caustics  was  shown  to  disappear  if  2D  Maslov  integrals  are  used.  Maslov  amplitudes  are  computed 
using  dynamic  ray  tracing  and  the  method  of  stationary  phase. 

The  theory  was  applied  to  study  Lg  propagation  in  southern  California.  Array  data  from  the 
NPE  at  the  Nevada  Test  Site  showed  that  Lg  disappears  at  epicentral  distances  between  250  km 
and  350  km.  By  using  ray  density  plots  it  was  shown  what  crustal  structure  is  needed  in  order  to 
show  this  behavior. 

The  algorithm  presented  here  concerns  acoustic  wave  propagation.  In  order  to  study  crustal 
waveforms  in  general,  and  not  only  Lg,  it  is  necessary  to  extend  the  algorithm  to  include  P-S  and 
S-P  conversion.  This  woidd  also  make  it  possible  to  study  the  change  of  P/Lg  amplitude  ratios. 
The  use  of  ray  density  plots  as  energy  plots  should  be  justified  using  synthetics.  An  excellent 
method  is  provided  by  the  Maslov  algorithm  presented  in  this  paper.  This  algorithm  also  makes 
it  possible  to  study  effects  of  3D  heterogeneity.  However  this  requires  the  existence  of  high  quality 
array  data. 
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Figure  2.  DeiiiJtious  of  the  unit  vectors  nnd  angle  used  in  the  kinematic  ray  equations 
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Figure  3.  Phase  space  plots  (left)  and  Lyapunov  exponents  (right)  for  9  rays. 
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all  sensors;  traveltlme  curves 


Figure  4.  NPE  data  from  the  Nevada  Test  Site.  The  reduction  velocity  is  6.9  km/s.  Traveltimes 
of  PmP,  SmS  and  midcnistal  multiples  are  shown.  Lg  arrives  after  SmS. 


Figure  5.  Source  and  station  locations  for  the  data  shown  in  figure  4.  Also  indicated  are  the 
major  tectonic  provinces  of  the  area. 
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